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SUMMARY 


To get a basic understanding of the physics of flowfields modeled by vortex 
filaments with finite vortical cores, systematic numerical studies of the 
interactions of two-dimensional vortices and pairs of coaxial axisymmetric circular 
vortex rings were made. Finite-difference solutions of the unsteady incompressible 
Navier-Stokes equations were carried out using vorticity and stream function as 
primary variables. Special emphasis was placed on the formulation of appropriate 
boundary conditions necessary for the calculations in a finite computational 
domain. Numerical results illustrate the interaction of vortex filaments, 
demonstrate when and how they merge with each other, and establish the region of 
validity for an asymptotic analysis. 


INTRODUCTION 

The vortex core structure, vortex merging, and the intersection of multiple 
vortex filaments have been of interest for many years. The flows are characterized 
by time-dependent and three-dimensional interactions. 

Mathematical models of vortex filaments have been frequently employed for the 
explanation of phenomena such as free jet flows [1] and aircraft wakes [2], 
Asymptotic analysis is well known to be valid if 6/£«l, where 5 is the effective 
viscous core size and £ is the characteristic flowfield length. When merging is 
about to occur, i.e. , as the ratio 6/ £ in the flow . becomes no longer small, the 
asymptotic theory may not be applicable. To continue the study of flowfields 
dominated by interacting vortex filaments, numerical solutions of the Navier-Stokes 
equations must be carried out. 

One of the main difficulties in numerical simulation of a flowfield in an 
unbounded domain is how to formulate appropriate boundary conditions necessary for 
solving the Navier-Stokes equations in a finite computational domain. In order to 
avoid a large computational domain, a method of specifying boundary conditions on a 
bounded domain with error control has been presented In [3]. Based on integral 
properties of the solution, this method allows the size of the computational domain 
to be reduced until the error introduced at the boundary approaches the same order of 
magnitude as that of the finite-difference approximation. Variations in the applica- 
tion of this efficient boundary condition method have been reported in [2,4-8]. 

To study the motion and merging of vortex filaments, the unsteady incompressible 
Navier-Stokes equations with aforementioned efficient boundary conditions are solved. 
The time variation of contour lines of constant vorticity illustrates the motion and 
merging of vortices. A much better picture of merging can be seen with graphs of 
trajectories of the locations of maximum vorticity and with graphs of the decay of 
maximum vorticity with respect to time. The practical limit for S/£ in the 
asymptotic theory is also shown. Cases included are (1) two, three, and four two- 
dimensional vortices; and (2) two coaxial axisymmetric circular vortex rings. 
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Subsequent sections of this paper contain a description of the governing 
equations, a brief discussion of the appropriate boundary conditions, an outline of 
the computational scheme, and numerical results demonstrating the characteristics of 
the vortical interactions. 
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SYMBOLS 

vector velocity potential 

effective radius (centered at the origin) of the vorticity distribution in 
two-dimensional cases or cross-sectional radius of a vortex ring 

error of approximate boundary data for A 

error of the finite-difference approximation 

error due to w = 0 on the boundary 

initial separation distance between two vortex rings 

size of the computational domain 

characteristic flowfield length 

number of two-dimensional interacting vortices 

cylindrical coordinates - 

Reynolds number, T/v 

initial centerline radius of a vortex ring 



time scale 

time-shift constant for an optimum single Lamb vortex 
initial time of numerical computation 
velocity vector 
inertial Cartesian coordinates 

rotating Cartesian coordinates fixed on one of the lines through the 
vorticity maxima » 

general position vector 

location of vortex maximum for asymptotic solutions 
circulation 

effective viscous core size 
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6 initial effective viscous core size 

o 

C scalar vorticity for Navier-Stoke solutions 

5 scalar vorticity for the sum of asymptotic solutions 

£* scalar vorticity for an optimum single Lamb vortex 

v kinematic viscosity 

p initial distance between the origin and an individual two-dimensional Lamb 

vortex 

£2 instantaneous angular speed of the x-axis 

a) vorticity vector 

u> o initial vorticity distribution 

GOVERNING EQUATIONS 

Incompressible laminar flowfields induced by initial vorticity distributions are 


governed by 

■ V*V = 0 (1) 

S t + (V • V) u> - (w*V)V = vAu (2) 

w = V x V (3) 

subject to the initial condition 

m(x,0) = “ Q (x) (4) 

and the boundary condition. 

V(x, t)+0 as |x|+» (5) 

For a flowfield dominated by vortices, the initia| data w are either of bounded 
support or decay exponentially in r, where r = |x|. As a result, the vorticity 

decays exponentially in r for all t > 0 , that is 

m(x,t) = 0 (e br ) ( 6 ) 


for some positive constant Equations (1) to (5) define an initial— value problem 

in an unbounded domain for to and V. 

The divergence-free condition of the velocity field, equation (1), allows v to 
relate to its vector velocity potential A by 

V = V x A (7) 

After imposing V»V = 0 , A then fulfills the vector Poisson equation 
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( 8 ) 


4 - 4 - 

AA = -a) 

The governing equations are now (2), (7), and (8) wit^h w and A as prime 

variables. At each instant, the vector velocity potential A, i.e^, the solution of 
equation (8) with boundary condition (5), can then be related to u> in an unbounded 
domain by the Poisson integral 
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A(x,t) 
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(9) 


Since numerical calculations can only tje performed in a bounded domain, 
appropriate boundary conditions for o> and A must be introduced. It is assumed 

that a) = 0 outside of the domain and that 

w = 0 on the boundary. (10) 


According to condition (6), the error of this approximation e is on the order of 
e - ^ with L as the size of the finite computational domain. The error of the 
finite-difference approximation e^ depends on the numerical scheme. By matching 
e^ with e , L can be chosen so that the approximate boundary conditon for w, 
equation (107, is acceptable. 


This is the approach employed for the two-dimensional problems in [9,10]. It is 
xjoted that for the two-dimensional or axisymmetric case the vector velocity potential 
A has only one component, i.e., the scalar stream function, and the vector vorticity 
w also has only one component, i.e., the scalar vorticity ?. The evaluation of a 
Poisson integral is extremely time Consuming even for two-dimensional problems. It 
will become worse for three-dimensional problems, because the number of operations 
is 0(N b ) for a grid with N points in each direction. Even if the Poisson 
integral is used to evaluate the boundary data for A and the Poisson solver [11] is 
applied to the interior domain, it is still very inefficient. The reason is that the 
number of operations for calculating the domain interior is only 0(N^log N) while 
the number for evaluating the boundary data is O(N^). The efficient boundary condi- 
tion specification technique used here was developed by Tine [3] so that the number 
of operations for establishing ^ the boundary data is 0(tP) and the error e^ 
introduced at the boundary for A is 0(L~ m ) with m>4. 


APPROXIMATE BOUNDARY CONDITIONS FOR VECTOR VELOCITY POTENTIAL 

To generate more accurate boundary data for A, the farfield behavior of the 
Poisson integral (9) is developed. When the vorticity decays exponentially in the 
farfield, A can be represented as a power series in r - * by expanding the 
denominator of the Poisson integral. The coefficients of the terms in the power 
series are moments of the vorticity distribution. The first moments and several 
linear combinations of higher moments are time-invariant. These results [3] are used 
to determine the boundary data for numerical solutions. Time-invariant moments of 
vorticity distributions can also be used to monitor the global accuracy of the 
finite-difference solutions. 

The error of approximate boundary data for A is e» = 0(L -m ) and m depends 
on the number of terms used in the power series. By matching e^ with e^, which is 
usually much larger than e , the optimal size of the computational domain is deter- 
mined. 
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The appropriate boundary conditions for A in the general three-dimensional case 
are presented in [3]. They can be reduced to the two-dimensional case [2, 3] or the 
axisymmetric case [3, 4] with proper modifications. 

COMPUTATIONAL SCHEME 

In two-dimensional cases, the effective radius a (centered at the origin) of the 
vorticity distribution is required to be sufficiently smaller than the size of the 
domain L, i.e., a « L, all the time. In axisymmetric cases, the vorticity distri- 
bution is considered to be concentrated almost entirely inside a torus whose cross- 

sectional radius a is much smaller than its centerline radius R, i.e., a « R. To 

study the merging of two vortex rings with very small core radii, the special case 
where a « L « R is considered. 

Since the vorticity is assumed to decay exponentially in distance, the initial 

data a) are generated by summing several individual Lamb line vortices (asymptotic 

solutions) and Lamb-type vortex rings for the two-dimensional and axisymmetric 
problems, respectively. The minimum distance between vortices is sufficiently larger 
than each individual core size so that asymptotic solutions are valid. 

The finite-difference scheme consists of the following steps. 

Vorticity integrals are evaluated using Simpson integration over the finite 
computational domain. 

The boundary data of A are calculated by applying the farfield approxima- 
tions metioned in the previous section. 

A fast Poisson solver is employed to determine A in the interior of the 
computational domain. 

The velocity field V is obtained by taking the curl of A. 

■4 

The vorticity field w is advanced in time by integrating the vorticity 
transport equation (2) in a conservative form using either the Duf ort-Frankel 
scheme or the ADI scheme. 

Steps (i) through (v) are repeatedly applied to advance the finite-difference 
solution in time. In order to attain a time-accurate transient solution to the 
unsteady flow problem, relatively small time steps must be taken. The time step 
should also be consistent with the truncation error due to spatial discretization <|nd 
with the error due to taking only finite terms of farfield approximations for A. 
This algorithm is accurate to the second order in space and time. 


(i) 

(ii) 

(iii) 

(iv) 

(v) 


RESULTS AND DISCUSSION 

The unsteady incompressible Navier-Stokes equations with efficient boundary 
conditions are solved for two-dimensional and axisymmetric interacting vortices. 
Examples are chpsen to demonstrate the merging of vortices and the practical limit of 
the asymptotic analysis [12-14]. 

For two-dimensional problems, the merging of M identical vortices located 

symmetrically on a circle of radius p is studied for M ** 2, 3, and 4. Under 
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inviscid theory, the vortices will move along the circle at a constant speed. Due to 
the decay of the vortical core and the enlargement of the effective individual core 
size, the M points of maximum vorticity will spiral inward faster and faster as the 
vortical cores overlap each other more and more. Eventually the M points of 
■maximum vorticity coincide at the origin and the constant vorticity lines begin 
circularizing around the origin. 


To carry out the Navier-Stokes solutions, the initial effective core size is set 
to be 0.5p for all cases. Initially the vorticity at the origin is of the order 
of e~^ and can be ignored. The initial core structure for each vortex is that of 
an isolated Lamb vortex created at t=0. The initial core size is chosen as the 
unit length scale and 6 /4v as the unit time scale; therefore, numerical compu- 
tation begins at t“t »1. 


From the Navier-Stokes solutions, the evolution of contour lines of constant 
vorticity for M vortices is presented in figure 1 for Re=T/v = 100 and M=3. 
Figures la-lh show the different stages of merging. At each instant, four contour 
lines at 0.995, e~*' 2 , e - *, and e ~^' 2 times its maximum value, £ (t), are 
shown. X and Y are the inertial Cartesian coordinates, while x and y are the 
coordinates rotating with one of the lines through the vorticity maxima. Figure la 
(displays the initial contour lines of the three nonoverlapping Lamb vortices. 

-3/2 

Figure lb shows that the outermost contour lines of C = e £ (t) for the three 
vortices are merging at t=1.60. Figures lc-le show the successive merging of the 
inner contour lines one by one at t=1.90, 3.11, and 4.02, respectively. For t>.62, 
the three maximum vorticity points merge at the origin and the character of three 
distinct vortices disappears. Figures If-lh then show the gradual circularization of 
the contours beginning from the inner ones as time increases. In figure lh, t=8.24, 
all four contour lines are nearly circles. The vorticity should then decay like a 
single Lamb vortex of strength MT with an appropriate age defined by an optimum 
asymptotic solution [10, 14]. 


When the initial vorticity distribution is defined by the sum of asymptotic 
solutions for nonoverlapping Lamb vortices, the solution 


C(x,t)= 


T_ 

4 TTVt 


M 

Z 

j = l 


-[ p 2 +x -2px cos( J ^~)2ir]/4vt 
M 
e 


( 11 ) 


is certainly quite accurate at the early stage. However, it is useful to find out 
how far the asymptotic theory can be extended into the overlapping stage by comparing 
the solution 4 with the ^Navier-Stokes solution 4 . From equation (11), the 
instantaneous angular speed of the x-axis (the first line through the maximum c) 
is obtained as 


r M -p 2 [ l-cos(ijp- 2 ir) ] / 2 vt 

&(t) -j Z {l-e M } 

4i j=2 


( 12 ) 


rw <v 

From equation (11), the location x(t) ^for ? (t) is found. Together with 

equation (12) the trajectory and decay of £ (t) are determined. 

infix 

Figure 2 illustrates that the difference between the decay of ? (t) obtained 

from the sum of asymptotic solutions and that of C (t) from Navier-Stokes 
solutions becomes noticeable only when merging is taking place. It is surprising 
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that the t curve follows very well with the c curve for M=4 even after 

, max max 

merging. 

Figure 3 indicates the trajectories of the locations of maximum vorticity for 
Navier-Stokes solutions (solid lines) and for the sum of M asymptotic solutions 
(dotted lines). The agreements of the trajectories between them become better as 
M increases; however, the trajectories of inviscid solutions (dash lines) remain as 
the circle of radius p and can never show the merging characteristics at all. 

It is found that even when the core size of an individual Lamb vortex Is on the 
order of the distance between vortices, the sum of asymptotic solutions is in good 
agreement with the Navier-Stokes solution. Also, the agreement remains good for a 
longer time with larger M. This can only be attributed to the fact that the non- 
linear convection terms become smaller and smaller as the number of vortices 
increases. Without the nonlinear terms, the sum of the asymptotic solutions is the 
exact solution. 

After merging, the Navier-Stokes solutions may be represented by an optimum 
single Lamb vortex defined by 


C (X,Y,t)= 


mt 

: k 

4w(t+t ) 


-(X 2 +Y 2 )/4v( t+t*) 


(13) 


where t* is a constant time-shift parameter obtained by letting the single vortex 
have the same polar moment as that of the initial data [10, 14]. This isolated 

vortex is created at t=-t*, i.e., when it is a point vortex. For the present initial 
vorticity distribution, it can be shown that 

* 2 

t * -jr^ + tg for all M. (14) 


It is equal to 5 when v=0.25. As shown in figure 2, the agreement between the 

decay of and that of £ becomes better and better as time increases, 
max max 

For axisymmetric problems, the cylindrical coordinates (R, Z) are used, where R 
is in the radial direction and Z is on the axis of symmetry. The coordinate frame 
translates in the Z-direction at an instantaneous inertial speed which is an average 
of the f lowfield velocity weighted by the local magnitude of the vorticity. 

Numerical examples for two initially identical coaxial vortex rings of radius 
R q =15 and of similarity vorticity distribution are carried out to show how and when 
they will merge into a single vortex ring. The initial core size 6 is used as the 
unit length scale. The unit time scale is chosen such that the initial time 
t =l/4v. The two rings are initially separated by a distance D Q . The Reynolds 
number based on the initial core size is defined as I*/ v where T denotes the 
strength of a vortex ring. 

Figure 4 illustrates a typical evolution of the vorticity for two right-moving 
interacting vortex rings with Re=754, v = 0.25, and D Q =6. Note that the displayed 
vorticity contours represent the vorticity distribution on a meridian plane cutting 
through the torus with only the contours centered at (R Q , Z) being shown. The center 
of the torus is under the contours. The first ring, moving ahead initially, is 
stretching its radius and slowing down. In contrast, the second ring, moving behind 
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initially, is contracting its radius and speeding up. The plane of the second ring 
passes over that of the first ring at t=1.83. At this instant, the two rings switch 
the roles of leading and lagging. It is observed that the two vortex rings cross 
over each other two more times at t=3.11 and 3.81. Finally, the two rings merge into 
"a single ring for t > 3.81. " p 

Figure 5 shows the decay of maximum vorticity for various initial separation 
distances D Q =3, 4 and 6 at a fixed Reynolds number of 251. As D Q decreases, the • 

vortical interactions become stronger so that the maximum vorticity decays faster and 
complete merging occurs sooner. Figure 6 displays the trajectories of the locations 
of maximum vorticity. The two rings pass over each other once, separate once, and 
then merge into a single ring in all three cases. 

Figure 7 depicts the decay of maximum vorticity for various Reynolds number 
Re=47, 189, and 754 at a fixed D Q of 6. As Re decreases, the viscous effects become 
larger so that the maximum vorticity decays faster and complete merging occurs 
earlier. It should be noted that the time is scaled by the Reynolds number. 
Figure 8 presents the effect of Reynolds number on the trajectories of the locations 
of maximum vorticity. The two vortex rings pass over each other three times before 
merging for Re=754, pass over each other only once before merging for Re=189, and 
merge into a single ring without passing over each other for Re=47. 

The difference of the maximum vorticities between the two rings decreases as Re 
or D q decreases. The reason is that the effect due to stretching of the ring 
diminishes as Re or D Q decreases. Consequently, both curves of the decay of the 
maximum vorticity for the two rings become closer to that of the corresponding two- 
dimensional results. 

CONCLUSIONS 

The unsteady incompressible Navier-Stokes equations have been solved for two- 
dimensional and axisymmetric interacting vortices. Numerical examples have been 
chosen to demonstrate the merging characteristics and establish the practical limit 
of the asymptotic analysis. 

For two-dimensional problems, it is found that even when the core size of an 
individual Lamb vortex is on the order of the distance between vortices, the decay of 
the maximum vorticity obtained from the sum of asymptotic solutions is in good agree- 
ment with that from the Navier-Stokes solutions and the agreement remains longer in 
time for larger numbers of vortices. When merging is about to take place, numerical 
solutions of the Navier-Stokes equations are required. After complete merging, the 
optimum single Lamb vortex solution can represent the Navier-Stokes solution better 
and better as time increases. The comparison of the trajectory of the locations of 
maximum vorticity for Navier-Stokes solutions with that for the sum of asymptotic 
solutions becomes better for larger number of vortices. On the other hand, inviscid 
theory can neither show the decay of the maximum vorticity nor the correct trajectory 
of the locations of maximum vorticity. 

For axisymmetric problems, numerical results show that as the Reynolds number or 
the initial separation distance decreases, the effect due to stretching of the vortex 
ring diminishes so that the maximum vorticity decays faster and complete merging 
occurs earlier. As the Reynolds number increases, the elapsed times for switching 
the roles of leading and lagging of the two vortex rings will increase. 
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